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Abstract 

We show that by a suitable change of variables, the derivatives of molecular 
integrals over Gaussian-type functions required for analytic energy derivatives can 
be evaluated with significantly less computational effort than current formulations. 
The reduction in effort increases with the order of differentiation. 

I. Introduction 

Analytic energy derivative methods have revolutionized the application of com- 
putational quantum chemistry to problems of chemical interest [1]. The location 
and characterization of stationary points on polyatomic molecular potential energy 
surfaces can be accomplished so much more efficiently using analytic derivatives 
than with techniques based on computing energies alone that the development and 
extension of analytic derivative methods has been one of the most active fields of 
methodological research in quantum chemistry in recent years. Given the gradi- 
ent and Hessian of the energy with respect to the nuclear coordinates, a variety of 
strategies have been developed that are guaranteed to converge to minima on po 
tential surfaces and that can efficiently locate other stationary points, particularly 
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transition states. These strategies can also be used to “walk” on surfaces from one 
minimum to another, thereby defining a reaction coordinate, and among the most 
elegant and conceptually illuminating studies of this sort are the investigations 
of Ruedenberg and co-workers on rearrangement reactions of small hydrocarbon 
species (see Refs. 2-5 and references therein). It is thus a great pleasure to dedicate 
this contribution to Professor Ruedenberg on the occasion of his 70th birthday. 

Of course, in order to perform such walks and optimizations it is imperative 
to evaluate the energy derivatives efficiently at the computational level of interest 
(Hartree-Fock or some correlated treatment). As noted above, much work has been 
performed in this area, and several reviews are available [1,6,7]. We shall concentrate 
here on a topic that ultimately affects the computational effort necessary to evaluate 
energy derivatives for any ab initio method that relies on a basis set expansion of 
Gaussian one-electron functions. 

Wave functions for polyatomic molecules are invariably expanded in a basis set 
that is centred on the various nuclei, and so in a calculation of the energy derivative 
of nth order with respect to the nuclear coordinates, up to nth-order derivatives of 
the one- and two-electron integrals are required. These derivative integrals can in- 
volve differentiation of the operators as well as differentiation of the basis functions, 
but the greatest computational problems arise from the differentiation of the basis 
functions. Like the evaluation of integrals over Gaussians [8,9], the calculation of 
integrals over differentiated Gaussians has been the subject of many investigations 
and numerous efficient computational schemes have been devised. In this work we 
show how the efficiency of derivative integral evaluation can be improved by some 
simple manipulations of variables. We shall briefly review the McMurchie-Davidson 
scheme [8] for computing Gaussian integrals and derivative integrals, and then show 
how a change of differentiation variables simplifies the formulas. 

II. Derivative Integral Formulas 

We shall expand the Gaussian charge distributions that appear m the integrals 
in Hermite functions, as described by McMurchie and Davidson [8] (see also Saun- 
ders [9]). Let us represent an unnormalized Cartesian Gaussian function centred 

at A by 

G iifc (f,a,A) = x^z'Xexpi-ar^), (!) 

where x A = x-A x ,e tc. We can consider one Cartesian direction, say x, represented 
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as 


( 2 ) 


Gi(x,a,A x ) = x l A exp{-ax A ). 

The overlap distribution of two such functions is expanded as 

a, 6, A x , B x ) = Gi(x,a,A x )Gj(x,b,B x ) 

= £ E t ij (a, 6, A x , £«) A*(x, P , P«), (3) 

t= 0 


where the Hermite function A t (x,p,-Pz) is defined by 



A t {x,p, p x ) = {d/dP x y exp(-px 2 P ) 

(4) 

with. 

P = -A+-B 
P V 

(5) 

and 

+ 

e 

II 

(6) 

The expansion 

coefficients E\ 3 (a,b, A x , B x ) are obtained from 



R X E\ 3 + {t + 1 )-^t+n 

2 p P 

(7) 

where 

Rx = B x 

(8) 


ana i 

E« a =exp(--Rl). (9) 

V 

Henceforth we shall not always list the arguments of the expansion coefficients or 
Hermite functions, but we wish to emphasize here that the expansion coefficients 
depend on a, b, and R x only, while the Hermite functions are independent of R x : 

i+i 

fty(z l a l M„B,) = £)E?'(a,&,fl*)A t (z,p,P,). ( 10 ) 

t=0 

In terms of the Hermite functions and expansion coefficients we can express a 
two-electron integral 



x\y k A z A exp(-ar^) x 3 B y l B z% exp {-br B ) 
x x l cyc z c' exp{- cr c) x Dyb z D exp(-dr 2 D ) dr 1 dr 2 


( 11 ) 
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as 


»+ j 


i'+3 


£ E V M,C„D.) 

t = 0 t'=0 

fc+i *'+»' , , 

x ^E^ l {a,b,A y ,By) Y, E u ' 1 (c,d,Cy,Dy) 


u=0 

m+n 


u '=0 
m' + n' 


X 


X 


y~^ E™ n (a,b,A z ,B z ) ^ E™ n (c,d,C z ,D z ) 


v=0 


v '=0 




( 12 ) 


where 


= JJ At(z,p,P*) E t '(x,q,Qx) &u(y,P,Py) E u '{y,<hQy) 
x A „(z,p, P z ) A V '( z ) 9? 


(13) 


and q and Q are defined analogously to p and P but for the second charge dis- 
tribution. Thus in practice we evaluate integrals over the Hermite function basis 
and combine those with the expansion coefficients to give integrals over primitive 
Gaussians. Some modifications to the form (12) are desirable from the point of 
view of efficiency, as discussed by Saunders [9], but for schematic purposes we can 
use (12). The first step, evaluation of the Hermite function integrals, is fast. The 
second step, which we can regard as a transformation from the Hermite function 
basis to the Cartesian Gaussian basis, is relatively time-consuming and is certainly 
more expensive than calculating the Hermite function integrals. Finally, if required, 
we combine these integrals with basis set contraction coefficients to give the final 
integrals. In fact, some of the expansion steps can be taken outside the contraction 

step, with a consequent improvement in efficiency. 

In a derivative integral we are interested in derivatives of dflij/dA x and 

dtlij/dBx for first derivatives, for example. Conventionally, we would differentiate 
the orbitals (2) first and then expand the overlap distributions of the differentiated 
orbitals analogously to fly above. For example, for the derivative with respect 

to A x we obtain 


Ollij __ ^ 

ft A ~ 

t= 0 


if A t . 


(14) 


4 



Note that the sum here is over more terms than appear in the undifferentiated charge 
distribution (3) — higher orders of differentiation would increase this summation 
range further. The new coefficients if are defined in terms of the coefficients if 

above by 

if = 2aE l t +1,j - iE\ ~ 1 ' J . ( 15 ) 

Analogous coefficients can be defined for higher orders of differentiation or for dif- 
ferentiation with respect to B x . In this approach, then, we compute derivative inte- 
grals using the same general scheme (12) as for undifferentiated integrals. Since the 
expansion of the differentiated charge distributions in Hermite functions is longer 
than for the undifferentiated distributions, the work required to transform from the 
Hermite function basis to the Cartesian Gaussian basis is greater. Further, as the 
order of differentiation increases this extra work becomes larger and larger. Hence 
this approach is not well-suited to higher derivatives. 

Let us instead consider differentiation with respect to the variables P x and R x , 


for which 

d 

dA x 

a d 
P dP x 

d 

+ dR x 

(16) 

and 

d 

b d 

d 

(17) 


dB x 

CD 

dR * 



We recall that the Hermite functions are independent of R x , while the expansion 
coefficients ate independent of J>„. Hence we can expect the expressions for the 
differentiated charge distributions to be simpler in terms of these variables, although 
we must eventually transform the derivatives back to the A„ representation. We 
obtain for the derivatives 


dfli j _ r ij dA t 

dP * 


E E't h-t+i 

t—0 


(18) 


ana • , • 

dfltj _ y- ^L-A ( 19 ) 

dR* h dR * 

Denoting dE\ j / dR x by E\ j '\ we obtain the expansion relation 

E i+1 ' jil = —E l t j; \ - -( R X E t ij;1 + E l t j ) + (t + 1 )E%\ (20) 

1 2p t_1 p 
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by differentiating the relation (7) above. 

We can make several important observations about these derivative formulas. 
First, the combination of expansion coefficients and Hermite functions m (18) above 
is over exactly the same range as the summation to give undifferentiated integrals: 
the only difference is that the degree of the Hermite function has increased by 
one. Hence the code required to evaluate this term is the same as required in the 
undifferentiated case, and the number of operations is also the same. (It is easy 
to see that this holds true in any order of differentiation for this term.) As we 
saw above, this is not the case if we differentiate with respect to the variables A x 
and B x , because then a linear combination of different degree Hermite functions 
and expansion coefficients appears. 

Second, calculation of the differentiated expansion coefficients E; 1 ’ 1 requires 
essentially the same code again as for the undifferentiated case, with the obvious 
addition of an extra term in the expansion relation, and a starting value 

eS ° !1 = -—R x E°°, ( 21 ) 

P 

obtained by differentiating (9). As we noted, the index range of the coefficients that 
are required is the same as that for the undifferentiated case, so the actual work 
required to combine Hermite function integrals and expansion coefficients does not 
increase. (The precomputation of the expansion coefficients themselves is of course 

a very rapid step.) 

Third, in the usual scheme the index range of the program loops over the 
variables t, u, v depends on the direction of differentiation (i.e., differentiation with 
respect to A x , A y , etc). Thus these loops must be executed with different ranges 
for each of the three directions for first derivative integrals, for example. With 
our transformation of variables, the loop index ranges become independent of the 
direction of differentiation, so the program logic is simplified and the overheads 
are reduced. We may also note here that this approach in no way diminishes the 
possibilities for vectorizing the calculation of the integral derivatives. Indeed, the 
simplifications to the program loop structure are likely to enhance these possibilities. 

Fourth, we can obtain an additional simplification as follows. Adding (16) 
and (17) we obtain 

?_ (22) 

dB x dP x dA x 
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Now, (in addition to saving one multiplication) this form of the expression for 
the derivative with respect to B x does not depend on the orbital exponents at 
all. Hence we can delay the transformation to the B x derivative until later in the 
calculation, for example, until after the contraction step, so that the time required 
for this variable transformation becomes negligible. This is most important for first 
derivatives, as in any order of differentiation only one term can be treated this way. 

In the case of higher derivatives there is a variety of terms to be considered but 
the scheme remains essentially the same. For example, the nth-order differentiated 
expansion coefficients with respect to R x are obtained from the recursion formula 


i ji n 


= —E 1 /' ? 

2 p 11 


f (R,E' t iin + nE? n -' ) + (t + mi? 
p 


(23) 


with starting values 


7 -i 00 ;nH~l 

^0 


P 


(24) 


and the identification 

£tj;° _ E ij (25) 

Higher derivatives of the Hermite functions with respect to P x (19) are trivially 
obtained. We note further that if the two charge distributions that appear in an 
integral are differentiated separately, the total savings is the product of the individ- 
ual reductions in work, since the two differentiations are independent. For multiple 
differentiation of the same charge distribution, we recall that by using our trans- 
formation of differentiation variables the summation range in the Hermite function 
to Cartesian Gaussian transformation is independent of the order of differentiation. 
Hence the savings increase as the order of differentiation increases, since in the 
conventional scheme the work required to accomplish this transformation increases 
substantially with the order of differentiation. In order to obtain an estimate of 
what savings are possible, we must also include an estimate of the effort required 
to transform back to the A X ,B X representation. We shall now present operation 
counts showing that it is always preferable to use our transformation of differenti- 

ation variables. 

In order to simplify the counting we consider only floating-point operations 
(multiplication and addition), which are weighted equally. In addition, m our count 
we have not taken advantage of the possibility of deferring transformation of some 
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derivatives until after contraction: in effect, we are counting operations only for 
primitive Gaussians and ignoring any additional savings that might accrue from 
moving manipulations outside the contraction step. If anything, neglecting this 
possibility favours the conventional approach to derivative integrals. 

We have listed operation counts for differentiation of SS, PP, and DD distri- 
butions in Table 1. We have not included the calculation of the Hermite function 
integrals, which is fast and contributes the same work to both cases, the conventional 
approach and our new scheme. Further, the transformation of the second charge 
distribution in the integral has also been excluded. We see that for the SS case the 
total operation count is not much affected by whether or not the transformation of 
variables is performed. However, for higher angular momentum functions there is 
a decided advantage to using the transformation of variables, and this advantage is 
clearly growing with the order of differentiation. As a further illustration of this, we 
note that for third derivatives of a PP distribution, for example, the conventional 
method would require 14 448 operations, while using the transformation of variables 
the work would be reduced to 8 340 operations: a savings of 42%. 

Finally, some other aspects of this scheme deserve comment. We note that 

d _ b d ~ a d (26) 

dR x p dA x p dB x 


Therefore, the operation ^ is not the same as the differentiation - £§ 7 • But 
if A and B coincide then the differentiation with respect to R x does not contribute 
to the energy derivative: only the differentiation with respect to P x contributes. 
This simplification is already used in the ABACUS program [10]. We also note 
that the use of translational invariance to reduce the computational labour is not 
affected by our transformation of variables: for first derivatives, for example, we 

have , T j T 

+ = 0, (27) 

dP x dQ x 

where I represents the two-electron integral in ( 11 ), from the use of translational 
invariance. 


Conclusions 

We have shown that by employing a transformation of differentiation variables, 
the work required to evaluate derivative integrals can be substantially reduced. The 
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advantages of our new approach increase both with the order of differentiation and 
with the angular momentum of the Gaussian functions involved. Savings will be 
obtained in the calculation of energy derivatives for any wave function that is ex- 
panded in a Gaussian basis. In particular, the economies obtained by applying these 
methods to the calculation of third or higher derivative integrals will be substantial. 
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Table 1. Operation counts for differentiation. 


SS PP 


First 

Derivatives 


Hermite/Cartesian Transformation 

12 

396 

P X ,R X to A X ,B X Transformation 

9 

81 

Total 

21 

477 

Conventional 

24 

672 


Second Derivatives 


Hermite/ Cartesian Transformation 
P X ,R X to A X ,B X Transformation 
Total 

Conventional 


42 

93 

135 


1 386 
837 

2 223 

3 678 


DD 


4 032 
324 
4 356 

6 144 


14 112 
3 328 
17 460 


150 


30 912 


